%% Stock Return Extrapolation, Option Prices, and Variance Risk Premium
% By Adem Atmaz
% This file generates Table 2: Option Prices
% It also requires 2 function files: RN_OptionPrices.m and RN_OptionPsi.m
% Default code generates Put prices for moneyness K=0.9*St
% For other moneyness levels change K to 1*St or 1.1*St in row 42
% For call option values change 'P' to 'C' in row 67
% Typically it takes less than 10 seconds to run

clear all;  clc;
format compact
format

% Parameter Values follow from  "Atmaz_Table_4_ParameterValues.m"
alpha=0.50;
theta=0.25;
mu=0.0212;
kappa=0.3567;
zeta=0.0044;
sigma=0.1625;
rho=0;
r=0.0056;

D0=100;
a=2.75;
vec_theta=[0, 0.25, 0.50, 0.75];
N_th=length(vec_theta);

% Implied Quantities
M_DELTA=((kappa+(rho*sigma)*(1+theta))^2)-((2+theta)*(1+theta)*sigma^2);
M_c=theta/(alpha*(1+theta));
M_b=-(2+theta)/(kappa+((rho*sigma)*(1+theta))+sqrt(M_DELTA));
M_sigma_1=(1+(rho*sigma*M_b))/(1-alpha*M_c);
M_sigma_2=(sqrt(1-rho^2)*sigma*M_b)/(1-alpha*M_c);
M_Var_con=(M_sigma_1^2)+(M_sigma_2^2);
M_Cov_con=(rho+(sigma*M_b))/(1-alpha*M_c);
M_corr_sv=M_Cov_con/sqrt(M_Var_con);

% Option contract quantities
St=100;             % Current stock price
Tt=3/12;            % Option Maturity
K=0.9*St;           % Option Moneyness. NOTE: vary this to 1 and 1.1 for other rows
V0=zeta/kappa;      % Long run mean of  fundamental  variance
vt=M_Var_con*V0;    % Current stock return variance is equal to its long run mean
cur_vol=sqrt(vt);   % Current stock volatility
vec_Xt=[-0.07, 0, 0.07];  % Varying vectors for columns Xt

OPrice=zeros(length(vec_theta),length(vec_Xt));
diff=zeros(length(vec_theta),1);
tic
for i=1:length(vec_theta)
    theta=vec_theta(i);
    DELTA=((kappa+(rho*sigma)*(1+theta))^2)-((2+theta)*(1+theta)*sigma^2);
    c=theta/(alpha*(1+theta));
    b=-(2+theta)/(kappa+((rho*sigma)*(1+theta))+sqrt(DELTA));
    sigma_1=(1+(rho*sigma*b))/(1-alpha*c);
    sigma_2=(sqrt(1-rho^2)*sigma*b)/(1-alpha*c);
    Var_con=(sigma_1^2)+(sigma_2^2);
    Cov_con=(rho+(sigma*b))/(1-alpha*c);
    kappa_v=kappa+(sigma*Cov_con);
    zeta_v=zeta*Var_con;
    eta_v=(theta/b)*(Var_con*(1-alpha*c));
    sigma_v=sigma*sqrt(Var_con);
    rho_sv=Cov_con/sqrt(Var_con);
    for j=1:length(vec_Xt)
        Xt=vec_Xt(j);
        OPrice(i,j) = RN_OptionPrices(kappa_v,zeta_v,eta_v,sigma_v,rho_sv,Tt,K,r,St,vt,Xt,alpha,'P');
    end
    diff(i)=100*((OPrice(i,1)/OPrice(i,end))-1);
end
toc
OPrice_diff=[OPrice, diff];
round(OPrice_diff,2)
 %% THE END

